Path-integral Monte Carlo method for the local Z 2 Berry phase 



o 

(N 
C 

(N 
(N 



i 

u 
-t— > 

CO 

c3 



i 

c 

o 
o 



> 



(N 



Yuichi Motoyama 1 and Synge Todo 2 

1 Department of Applied Physics, University of Tokyo, Tokyo 113-8656, Japan 
2 Institute for Solid State Physics, University of Tokyo, 
7-1-26-R501 Port Island South, Kobe 650-0047, Japan 
(Dated: January 23, 2013) 

We present a new loop cluster algorithm Monte Carlo method for calculating the local Z2 Berry 
phase of the quantum spin models. The Berry connection, which is given as the inner product of 
two ground states with different local twist angles, is expressed as a Monte Carlo average on the 
worldlines with fixed spin configurations at the imaginary-time boundaries. The "complex weight 
problem" caused by the local twist is solved by adopting the meron cluster algorithm. We present the 
results of simulation on the antiferromagnetic Heisenberg model on an out-phase bond-alternating 
ladder to demonstrate that our method successfully detects the change in the valence bond pattern 
at the quantum phase transition point. We also propose that the gauge-fixed local Berry connection 
can be an effective tool to estimate the quantum critical point precisely. 
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In low-dimensional systems or frustrated magnets, 
strong quantum fluctuations often disrupt the classical 
long-range orders. For example, the ground state of the 
5=1 antiferromagnetic Heisenberg (AFH) chain is not 
the classical Neel state but the Haldane state which 
is non-magnetic (5 to t = 0) and gapful. This state is 
understood in terms of the valence bond solid (VBS) 
picture [2]. The ground state of 5 = 1/2 AFH lad- 
der is another spin-gap state - the spin dimer state. 
While these models exhibit quantum phase transitions 
between a spin-gap state and another, these states are 
purely quantum state and thus they cannot be charac- 
terized by any classical order parameters. In decades, 
such states and phase transitions have been studied by 
using topological order parameters. In spin systems,the 
string order parameter Q and its generalizations 0-11] 
have been used to describe the Haldane and the dimer 
phases. Although the string order parameters have been 
used widely especially in numerical studies, they work 
only in one-dimensional systems. 

The local Z2 Berry phase, recently proposed by Hat- 
sugai 0, is one of the topological order parameters. This 
is a phase that the ground state wave function acquires 
under the twist of a bond on the lattice. The value of 
the local Z2 Berry phase is defined on each bond, pla- 
quette, and so on, and quantized to or tt (mod 27r) 
when the system has a finite gap and the antiunitary 
symmetry such as the time reversal symmetry. This or- 
der parameter is stable against any perturbation unless 
the gap collapses. Therefore, one can catch the quantum 
critical points by seeing the change of the spatial pat- 
tern of the Berry phases even at finite-size simulations. 
In the AFH models, for example, the value of the local 
Berry phase on a bond is 7r or depending on whether a 
singlet spin pair exists on the bond or not, respectively. 
Other models that have studied by this order parameter 
arc itinerant electron systems Q and Kondo singlets [t| . 



It is considered that this topological order parameter can 
be applied to two- or higher-dimensional systems as well 
as one-dimensional ones. The Berry phase and the local 
Z 2 Berry phase are related [? ? ] to the topological 
entanglement entropy, which is another tool to detect 
nontrivial topological phase [? ? ] . 

The local Z2 Berry phase, however, has been calculated 
numerically only by the exact diagonalization and thus 
the simulation size was strongly restricted (in S = 1/2 
AFH models, up to about 40 spins). To verify the con- 
jecture, that is, the local Z2 Berry phase also works for 
higher-dimensional systems, and to apply the order pa- 
rameter for attacking non-trivial quantum phase transi- 
tions, an unbiased numerical method, such as the path- 
integral quantum Monte Carlo (PIQMC) [10], that can 
treat larger systems is essential. 

The PIQMC calculation for the Berry phase has never 
been done since this has two obstacles. The first obstacle 
is that the Berry connection is defined as the inner prod- 
uct of wave functions of two ground states which cannot 
be accessed directly by the ordinary PIQMC. Second, a 
complex- valued weight appears due to the local twist in 
the interaction. This causes the "complex weight prob- 
lem," which is as difficult to deal with as the notorious 
negative sign problem. In the present paper, we propose 
a new procedure to calculate the local Z2 Berry phase 
by the PIQMC with the loop algorithm [ll|, Qj], and 
present some results of PIQMC simulation for 5=1/2 
AFH model on out-phase bond-alternating ladder. Our 
method works for any models in arbitrary dimensions as 
long as the original (untwisted) model has no negative 
sign problem. Extension to other models is straightfor- 
ward. We also propose a possibility that the gauge-fixed 
local Berry connection can be used as another topological 
order parameter. 

First, we define the local Berry phase of AFH models. 
To define the local Berry phase at bond (k£), one twists 
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this bond, that is, replace the term in the Hamiltonian 

Sk s l+ S k s i b y ei6 st s i+ e ~ ies k s t' where S± is the 

spin ladder operators and 9 is a twist angle. The local 
Berry connection is calculated as A{9) = {ip(9)\dg\4>(9)) 
by using the ground state \ip(8)), and the Berry phase is 
then defined as 
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To calculate this numerically, the integral needs to be dis- 
cretized as 7 = liniM-^oo J^jL^ 1 ar S W^IV^j+i)), 
where 9j = 2%j /M . 

The first obstacle in the PIQMC is how we calcu- 
late these inner products. As explained below, we fur- 
ther deform the inner product Bj = arg (ip(9j)\ip(9j + i)) 
into a Monte Carlo (MC) expectation form: Bj = 
ar §Ss ( ^( s )^ / ( s )/Ss^( s )' wriere s is a configuration 
index, O(s) is some c-number observable and W(s) is 
some positive weight function. The ground states can 
be represented by the projection method as \tp(9)} = 
lim^oo Cexp [— (3H(9)] \<fi(9)), where C is a real nor- 
malization factor, \4>(9)) is an arbitrary state which is 
nonorthogonal to the ground state, and (3 is a projec- 
tion parameter. Although the ^-dependence of \4>{9)) 
can be chosen arbitrary as long as |</>(0)) = |0(27r)) ; in 
the present paper we use the classical Neel state \4>) for 
all 9 for simplicity. The normalization factor C can be 
neglected since we need not the absolute value but the 
argument. Thus, the projection method leads us to 
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As well as the conventional continuous-time PIQMC [12| , 
we expand the exponential operators by path integral and 
introduce the worldline representation as an MC configu- 
ration and a complex valued weight function [see FIG. Q] 
(a)], by which Bj can be represented as arg^ s Wj(s). Fi- 
nally, reforming the complex weight to an amplitude and 
a phase factor, and dividing by the (real) normalization 
factor \ w \ lead us to the following MC form: 

w j (f) 

(3) 



Bj = arg ■ 



\Wj(s)\ 
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E s K'0)l 

Note that this deformation is justified since it does not 
change the argument. 

There are two differences from the conventional 
PIQMC: First, while the conventional PIQMC has the 
periodic imaginary-time boundary condition due to the 
trace operation, our method has the fixed one due to the 
gauge fixing [14[. Second, in our simulation the Hamil- 
tonian depends on imaginary time, H(t) = H(9j+i) or 
V.{9j), if t < (3/2 or r > (3/2, respectively. Since the 
effect of twist appears only as the phase of off-diagonal 
elements of the Hamiltonian, the phase factor is simply 

written as 
1 u ■ 

'9j (n£ - «u ) + Mj+i («l - "l 
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(a) spin configuration S 



(b) loop configuration (j 
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FIG. 1. (a) One of the typical worldline configurations of 
the S = 1/2 AFH model. Solid (dashed) vertical lines and 
upward (downward) arrows denote S z — 1/2 (—1/2) states. 
When one twists (hi) bond, = — 1, riy = = 
0. (b) One of the loop configurations compatible with the 
spin configuration (a). Current n±(c) of loop c (thick line) is 
n± (c) = (1 — 0) ± (0 — 0) = 1 and so this is the meron cluster 
when the twist angle 9 = ir. 



where (n^f ) is the number of off-diagonal operators 
S^S^ on the twisted bond (k£) at imaginary time r < 
(3/2 (t > (3/2), and (n^) is the number of off- 
diagonal operators S^S^. Thus, the MC estimator of 
the Berry phase is obtained as 



7 = ^arg( exp 



iQj (n^ - n^) + i9 J+1 (nj? - ri£) 



(5) 

where (• ■ ■ )o denotes the MC expectation value with re- 
spect to untwisted system. Furthermore, we can take 
the continuous limit M — > 00 by introducing variables 
9'j = (0j + 9 j+ i)/2 and A0 = 9 j+x - 9 3 = 2tx/M: 
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where n± = (n£^ — n \^) i ("u — n \j) an d we used the 
fact that (e l9n+ ) is real due to imaginary-time rever- 
sal symmetry. By comparing it with the definition of 
the Berry phase [eq.(JTJ|], the Berry connection with the 
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present gauge can be calculated by 



A(6) 
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As the projection parameter (3 becomes larger, the 
width of the distribution of n+ gets broader and so both 
(e m+e ) and (n_e m+e ) will be exponentially smaller for 
8 7^ 0, as e~ a P d with some constant a. This means 
that the relative error of the ratio of these two quanti- 
ties [eq. ((Tj)] will also diverge exponentially. Especially, 
at 8 = 7r this is nothing but the negative sign problem, 
which wc commonly encounter for the frustrated AFH 
model. The difference from the ordinary fully frustrated 
models (e.g. triangular AFH, kagome AFH) is that the 
local twist does not introduce the frustration in the di- 
agonal terms but only in the off-diagonal ones, and this 
fact fortunately enables us to apply the meron cluster 
algorithm 15], which is used to solve the negative sign 
problem for some fermion systems. Especially, the meron 
cluster algorithm completely solves the sign problem for 
the most severe case, that is, 8 = ir as seen below. 

In the following we will construct the improved esti- 
mators 16| for the denominator and the numerator of 
the Berry connection [eq. ([7])]. and apply the meron 
algorithm to this problem. The loop algorithm up- 
dates worldline configurations via loop configurations 
(see FIG.[T|). We can calculate observables in terms of the 
loop configuration by tracing out spin variables, which is 
referred to as the "improved estimator." First, we define 
local currents ro+(c) and n_(c) of each loop c as 



n± 



(c)=(nl(c)-ni(cj) 



± 



to 
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where n£ (nfj is the number of downward-upward 
(upward-downward) turns of the loop on the twisted 
bond at r < (3/2 and n\j {n\j) is the one at r > (3/2. 
By using n±{c), we can express the expectation values of 
the denominator and the numerator in eq. ([7]) as 



(e m+ x= n 



1 + e i6n +( c ) 



(9) 



and 




respectively, where (• • • ) g means the average over the 
loop configurations, and summation and production are 
taken over closed loops. When there are loops with 
8n+(c) — 7r mod 2ir, the contribution to the denomina- 
tor of eq. ([7]) vanishes. Such loops are called as "meron 
clusters." For example, the loop denoted by thick line in 
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FIG. 2. Out-phase bond-alternating ladder (0 < <5 < 1) . 
Singlet pairs are on strong leg bonds (thick lines) or on rung 
bonds (dashed lines) if J' < J' c {8) or J' > J' c (8), respectively. 



FIG. QJb) nas n +i c ) = n_(c) = 1, and thus this loop is 
a meron cluster if 6 = tt. When there are more than one 
meron cluster the numerator of eq. ([7]) is also zero. In 
order to sample loop configurations with non-zero contri- 
bution more frequently and reduce the statistical error of 
the Berry connection, we adopt the histogram reweight- 
ing with Wang-Landau optimization .17] with respect to 
the number of the meron clusters. Since — iA{8) is an 
even function due to the imaginary-time reversal sym- 
metry and has a period 2tt, we fit the numerical data for 
the discrete values of 8 by Y^k a k cos(k8) and obtain the 
final estimate of the Berry phase as 2nao by integrating 
the fitting curve. 

To demonstrate the present method we calculate the 
local Berry phase of S = 1/2 AFH model on out-phase 
bond-alternating ladder (FIG. [2|. The Hamiltonian of 
the model is 

L 

(ii) 

+j(l-(-l)^) S,. r S,.,.i 
+J'Si.j ■ S 2 ,j], 

where L is the ladder length (and so the number of sites is 
N = 2L) and S^j stands for the 5 = 1/2 spin operator 
on the j-th site of the «-th leg. The boundary condi- 
tion along the ladder is periodic, that is, S^l+i = S%i 
(i = 1,2). This model exhibits a quantum phase tran- 
sition at J' = J' c {8), where the pattern of singlet pairs 
changes globally; valence bonds are on strong le g bo nds 
or on rungs for J' < J' c or J' > J' c , respectively [18(. In 
the present study, we fix J = 1 and S — 0.5, for which 
the quantum critical point has been estimated as J' c ~ 
1.2 [19(. We choose M = 32 and generalized ensemble 
weight for loop configuration g as w(g,n m ) = w(g)/e nm 
for 8 ^ 7r while wc use Wang-Landau method for 8 = tt, 
where n m is the number of meron clusters. 

Figure shows the result for system sizes L — 8, 16, 32 
and projection parameter j3 = 2L. As J'/ J becomes 
larger, the local Berry phase on the strong leg bonds 
decreases from tt to while the Berry phase on the rung 
bonds increases from to tt. The projection parameter 
is not large enough to obtain ground states at near the 
critical point and thus these curves are not step functions. 
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FIG. 3. (color online) The local Berry phase on the out-phase 
bond-alternating ladder with system size L = 8 (squares), 16 
(circles), and 32 (triangles) on the leg bonds (closed red sym- 
bols) and the rung bonds (open blue symbols). The projection 
parameter j3 is 2L. 
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Finally, we point out numerical evidence that the 
gauge-fixed local Berry connection itself can be a topolog- 
ical order parameter. In fact, the Berry phase curves for 
several system sizes cross at near the real critical point, 
and thus there is no size-dependence of the Berry phase 
by Monte Carlo calculation at the critical point. We ex- 
pect that the gauge-fixed local Z 2 Berry connection does 
not also depend on system size at the critical point for 
all twist angle 9. Note that under the gauge transforma- 
tion, \ip(9)) —> e tx ^\ , ip(9)), the Berry connection varies 
as A{9) — > A(9) + idex(9), where \ 1S some arbitrary 
periodic real function such as x(6>) = x(@ + 2tt). This 
means that the J'- A curve with fixed 9 only shifts verti- 
cally by a constant under the gauge transformation and 
the crossing point J' c does not change. Thus, it is ex- 
pected that the gauge-fixed Berry connection curves will 
be on one universal function / by FSS transformation, 
A(J',N) = /((J' - J' C )N 1 ' ,J ). Figure [5] shows the re- 
sult of the Bayesian FSS analysis [20( of the imaginary 
part of the local Berry connection at 9 = 0. The system 
sizes are L = 128, 192, 256 and the projection parameter 
is (3 — 2L. The scaling parameters are J' c = 1.2268(1) 
and v = 0.738(8) for the leg twist and J' c = 1.2265(1) 
and v = 0.747(10) for the rung twist. These estimates 
for the critical point agree with the result of the FSS of 
staggered susceptibility. The results of several methods 
are summarized in TABLE [J] 

Calculating the Berry phase for one parameter requires 
some independent simulations. In order to suppress 
complex weight problem one has to use the reweight- 
ing method for the number of meron clusters, which is 
a global variable and so the simulation is heavy and dif- 
ficult to be parallelized. The local Berry connection for 
9 = 0, on the other hand, can be calculated from one 



FIG. 4. (color online) The estimation of critical point of out- 
phase bond-alternating ladders obtained by the local Z2 Berry 
phase on leg bond (circles), rung bond (triangles), and their 
crossing point (squares). The horizontal line, J' c /J = 1.2268 
is the FSS result of staggered susceptibility. 



In this case, however, since the energy gap remains finite 
except at the critical point, the curves convergence to 
step functions as L and /? become larger. 

We define three points to estimate critical points: 
f c lcs (L) and j' c mns {L) are the points where 7 lc s = tt/2 
and 7 run s = tt/2, respectively, and J c CTOSS (L) is the one 
where 7 lcg = <-y run s. The critical point in the thermo- 
dynamics limit, J' cl is estimated by size extrapolation of 
J'JL) for lattice sizes up to L = 32; = 1.2281(18), 
j c run g = 1.2282(18), and j' c cmss = 1.2266(6) (Fig. |4}. 
These results are consistent within statistical errors with 
the independent finite-size scaling (FSS) analysis for the 
staggered susceptibility, J' c / J = 1.2268(2). 
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FIG. 5. (color online) The FSS plot for the imaginary 
part of the gauge-fixed Berry connection, with system size 
L = 128, 192, 256. The downward to the right curve and the 
upward one are the leg twist and the rung twist Berry con- 
nection, respectively. 
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TABLE I. The quantum critical point J' c / J of out-phase 
bond-alternating ladder with dimerization 8 — 0.5 obtained 
by each method. 



method 


^max 


J'JJ 


Berry phase (leg twist) 


32 


1.2281(18) 


Berry phase (rung twist) 


32 


1.2282(18) 


Berry phase (cross) 


32 


1.2266(6) 


Berry connection (leg twist) 


256 


1.2268(1) 


Berry connection (rung twist) 


256 


1.2265(1) 


FSS of susceptibility 


256 


1.2268(2) 



simulation without the meron cluster algorithm and the 
reweighting method, and thus it is easy to calculate it 
and at the same time simulation can be parallelized. 

In conclusion, we presented the procedure to calculate 
the local Z2 Berry phase by the PIQMC and also pro- 
posed the gauge-fixed local Berry connection can be used 
as another topological order parameter. We calculated 
these order parameters of the AFH ladder models for a 
demonstration, and the result for the thermodynamics 
limit are consistent with the result by the conventional 
methods. 
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